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Sub-Grid  Models  for  Large  Eddy  Simulation  of  Turbulent 
Combustion 


1.1  Introduction 

Large  eddy  simulation  (LES)  increasingly  is  being  used  to  study  turbulent  reactive  flows  [1, 
2,  3,  4].  While  LES  is  computationally  more  expensive  than  ensemble  averaged  approaches, 
the  resources  that  LES  requires  are  applied  in  a  particularly  efficient  manner.  These  re¬ 
sources  are  used  to  resolve  the  large  energetic  scales  of  turbulence.  Large  scale  structures 
in  general  cannot  be  expected  to  remain  statistically  homogeneous,  and  ensemble  averaged 
models,  therefore,  have  difficulty  describing  them.  In  reactive  flows  in  particular,  large  scale 
inhomogeneities  are  very  important  because  they  influence  the  mixing  processes  that  affect 
the  evolution  of  combustion. 

The  inherent  resolution  of  the  LES  technique  reduces  the  amount  of  work,  relative  to 
ensemble  averaged  techniques,  that  computational  sub-grid  models  must  perform.  But 
the  need  for  accurate  sub-grid  scale  flow  descriptions  has  not  been  eliminated.  Rather, 
models  that  can  account  for  the  behavior  of  the  smallest,  unresolved  turbulent  scales  remain 
fundamentally  important.  In  spite  of  the  significant  progress  that  has  been  made  in  the  field 
of  LES,  many  questions  regarding  the  large  eddy  simulation  of  reactive  flows  in  particular 
still  exist.  Furthermore,  although  many  sub-grid  approaches  for  reacting  flows  have  been 
proposed,  universally  accurate  combustion  models  continue  to  remain  elusive.  Part  I  of  this 
research  project,  therefore,  sets  the  goal  of  analyzing  and  developing  improvements  to  a  set 
of  the  most  promising  sub-grid  models  for  LES  of  turbulent  combustion.  These  models  are 
based  primarily  on  asymptotic  flamelet  ideas  that  use  single  dimension  chemistry  solutions 
to  describe  multi-dimensional  turbulent  processes. 

Although  this  mapping  is  difficult  to  perform  accurately  in  itself,  an  additional  difficulty 
is  introduced  when  these  asymptotic  models  have  to  be  applied  simultaneously  to  multi¬ 
regime  combustors.  In  the  following  sections  work  will  be  presented  both  on  individual 
asymptotic  mapping  methods  and  on  how  to  combine  these  regime-specific  methods  appro¬ 
priately  in  realistic  device  simulations.  This  work  will  be  introduced  through  a  brief  review 
of  the  proposed  project  objectives. 

1.1.1  Review  of  Objectives 

Modern  premixed  turbulent  combustion  models  are  known  to  be  subject  to  a  number  of 
deficiencies.  For  example,  it  is  especially  difficult  to  describe  the  propagation  speed  of 
turbulent  premixed  flames  correctly,  because  flame  propagation  is  highly  dependent  on  what 
happens  at  the  smallest,  unresolved  flow  scales  [5,  6].  Second,  efforts  to  model  underresolved 
premixed  flames  such  as  the  level  set  approach  often  lead  to  inappropriate  descriptions 
of  premixed  flame  structure.  These  efforts  consequently  lead  to  incorrect  descriptions  of 
density  being  fed  to  a  flow  solver.  The  objective  of  the  LES  combustion  modeling  portion 
of  this  work  was,  therefore,  to  improve  flame  speed  models  and  to  lay  the  foundation  for  a 
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new  flame  structure  model. 

Work  in  the  area  of  non-premixed  turbulent  combustion  focused  on  improving  models 
for  the  sub- grid  variance  of  scalars.  The  variance  of  the  mixture  fraction  scalar,  in  particu¬ 
lar,  is  important  in  non-premixed  combustion  because  the  higher  order  moments  of  sub)- grid 
mixture  fraction  distributions  play  a  leading  role  in  determining  how  combustion  evolves. 
While  both  transport  equation  and  dynamic  algebraic  variance  models  exist,  this  project 
focused  on  improving  the  dynamic  procedure  relevant  for  algebraic  models.  Algebraic  mod¬ 
els  are  used  very  widely,  and  improvements  in  this  area,  in  turn,  may  help  improve  the 
modeling  of  unclosed  terms  in  transport  equation  approaches. 

Premixed  and  non-premixed  behavior  often  both  are  observed  in  modern  combustor 
designs  [1,  7].  If  the  advantages  of  flamelet  models  are  to  be  retained  in  these  multi-regime 
environments,  methods  of  selecting  whether  to  map  premixed  or  non-premixed  chemistry 
solutions  into  a  flow  solver  must  be  developed.  Until  recently,  the  flame  index  was  the  only 
available  selection  method.  The  flame  index  suffers  from  a  number  of  drawbacks,  however, 
such  as  its  inattention  to  information  regarding  chemical  time  scales.  Another  objective  of 
this  work  was,  therefore,  to  develop  an  improved  combustion  regime  index  that  would  allow 
flamelet  methods  to  be  combined  in  multi-regime  flows. 

1.2  Premixed  Turbulent  Combustion 

In  typical  large  eddy  simulations  of  premixed  combustion,  flame  fronts  exist  on  sub-grid 
length  scales.  Consequently,  reactive  scalar  variables  cannot  be  transported  accurately. 
Specifically,  the  changes  in  density  that  occur  at  the  flame  front  are  resolved  on  only  one  or 
two  mesh  cells.  These  poorly  resolved  density  changes  lead  to  large  numerical  errors  that 
affect  the  solution  of  all  reactive  scalar  transport  equations.  The  level  set,  or  G-equation, 
premixed  combustion  model  was  proposed  in  response  to  this  problem  [6,  8].  This  model, 
rather  than  resolving  changes  through  the  flame  front,  tracks  the  level  set  surface  that 
describes  the  front. 

The  accuracy  of  level  set  based  models  is  known  to  be  particularly  dependent  on  the 
numerical  methods  that  are  used.  In  response,  a  study  of  the  accuracy  of  numerical  level  set 
methods  in  the  context  of  premixed  LES  was  conducted.  A  number  of  canonical  level  set  test 
cases  exist,  and  many  of  these  have  been  used  to  validate  the  combustion  models  that  are 
discussed  below.  The  results  of  such  tests,  however,  cannot  be  used  readily  to  understand 
what  happens  in  turbulent  flow  fields  where  front  propagation  is  a  multi-scale  phenomenon. 
To  overcome  this  problem,  an  adaptive  mesh  refinement  technique,  the  Refined  Level  Set 
Grid  (RLSG)  approach,  was  implemented  and  applied  to  a  turbulent  flame. 

In  this  technique,  a  secondary  computational  mesh  is  created  locally  around  the  flame 
front  of  interest.  This  mesh  is  resolved  more  highly  than  the  baseline  flow  solver  mesh 
and  adapts  to  the  flame  front  as  it  moves.  The  highly  resolved  mesh  provides  a  means  of 
transporting  a  significantly  more  accurate  level  set  solution  than  the  solution  from  the  flow 
solver  mesh.  The  influence  of  numerical  error  on  the  flame  location  can  then  be  assessed 
by  running  with  increasingly  refined  level  set  meshes. 

If  the  numerical  errors  that  occur  when  the  level  set  is  solved  on  the  flow  solver  mesh  are 
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significant,  then  the  more  accurate,  refined  flame  solutions  should  display  different  charac¬ 
teristics.  Additionally,  the  refined  solutions  should  converge  to  the  correct  solution  as  the 
refined  mesh  resolution  is  increased.  Application  of  the  RLSG  to  a  turbulent  bunsen  flame, 
however,  showed  that  the  flame  front  solution  remained  roughly  constant  as  the  resolution 
increased.  Moreover,  the  solution  differences  that  did  appear  were  small  in  relation  to  the 
differences  that  developed  as  a  function  of  heat  release  modeling  choices.  For  example,  an 
adjustment  of  the  variance  used  in  the  probability  density  function  (PDF)  describing  the 
likelihood  of  finding  burned  gas  around  the  front  impacted  the  flame  location  more  than  the 
level  set  mesh  refinement  did.  The  RLSG  study,  therefore,  demonstrated  that  a  level  set 
framework  was  sufficiently  numerically  accurate  to  model  turbulent  premixed  flames.  While 
numerical  errors  were  created,  they  were  not  the  most  dominant  source  of  error.  Attention 
was  given,  therefore,  to  modeling  the  turbulent  burning  velocity,  premixed  heat  release,  and 
flame  structure. 

1.2.1  Dynamic  Turbulent  Burning  Velocity  Model 

Because  of  the  poor  flame  resolution  that  is  inherent  in  premixed  LES,  the  propagtion  speed 
of  a  flame  front  cannot  even  be  captured  moderately  well  by  a  filtered  transported  equation. 
This  speed,  therefore,  has  to  be  prescribed.  A  large  variety  of  models  for  this  propagation 
speed  have  been  proposed,  but  most  depend  on  a  series  of  empirically  defined  coefficients  [8, 
9,  10 j .  In  an  effort  to  improve  turbulent  flame  speed  predictions,  a  dynamic  procedure  was 
developed  for  calculating  these  coefficients  [5].  This  procedure  was  formulated  particularly 
to  be  consistent  with  the  level  set  approach  but  could  be  modified  easily  to  fit  other  models. 
By  analogy  to  the  dynamic  Smagorinsky  viscosity  model,  it  was  thought  that  a  dynamic 
procedure  that  accounts  for  local  information  could  describe  how  small  scales  influence  front 
propagation  more  accurately.  Another  advantage  of  the  dynamic  procedure  is  that  it  limits 
the  ability  of  model  developers  to  affect  flame  behavior  artificially  through  constant  tuning. 

The  derivation  of  the  dynamic  turbulent  burning  velocity  model  begins  by  formulating 
an  equation  describing  the  evolution  of  a  flame  front.  When  a  flame  front  is  defined  as 
an  isocontour  of  a  progress  variable  ( C  =  Co),  the  appropriate  equation  is  produced  by 
multiplying  the  progress  variable  transport  equation  with  a  delta  function, 


(i) 


This  delta  function  can  be  written  equivalently  as  the  derivative  of  a  heaviside  function, 
8(C  —  Co)VC  =  V[H  ( C  —  Co)]  and  can  be  combined  with  the  other  terms  in  the  transport 
equation.  Some  manipulation  produces  an  expression  describing  the  transport  of  a  heaviside 
function  located  at  the  flame  front, 


d[H(C-Co)\  ,  d[H  (C  -  C0)] 
di  +Uj  dxj 


=  Dk |V  [H  ( C  -  Co)] I  +  sL,c0|V  [H  (C  -  Co)] I,  (2) 


where  the  propagation  speed  of  the  Co  surface,  sl,C0>  has  been  introduced  to  account  for 
the  source  terms.  Unlike  level  set  formulations,  this  equation  is  valid  everywhere  in  space 
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and  not  just  at  the  flame  front  surface.  This  property  is  important  because  it  means  that 
a  standard  LES  volumetric  filter  can  be  applied  to  the  equation.  Making  the  notationally 
convenient  substitution  Q  =  H  (C  —  Co)  and  then  filtering  gives 


f +  ”^  =  (Dk)t  |ve|+JT,c.|vs|, 

(3) 

where 

«T, Co  |V0|  =  SL.ColV^I, 

(4) 

and 

(m)T  |vs|  =  dk\vg\, 

(5) 

represent  the  introduction  of  definitions  for  the  filtered  burning  velocity,  ^t,Cq ,  and  the 
filtered  propagation  speed  due  to  curvature,  (Dk)t. 

Equations  (3)-(5)  are  significant  because  they  can  be  used,  after  some  further  manip¬ 
ulation,  to  determine  the  form  of  the  terms  that  the  turbulent  burning  velocity  describes 
explicitly.  Additionally,  when  both  filtering  and  test  filtering  are  used,  they  provide  the 
framework  for  producing  a  dynamic  equation  that  can  be  solved  to  determine  the  value 
of  turbulent  burning  velocity  model  coefficients.  For  example,  since  filters  commute  with 
temporal  derivatives,  the  application  of  a  broad  test  filter  to  Eq.  (3)  produces  a  transport 

equation  for  Q  (the  arrow  operator  denotes  the  test  filter).  If  a  filter  and  a  test  filter  instead 

are  applied  to  Eq,  (2)  in  a  single  operation,  a  different  form  of  the  Q  equation  is  realized. 
Because  these  two  equations  appear  in  different  forms,  they  can  be  subtracted  from  one 
another  to  determine  an  expression  relating  turbulent  burning  velocity  models  associated 
with  different  filter  widths.  This  expression  can  be  written, 

+  |vg|  =  ((^)^  +  ^=pV,u)  |V?|,  (6) 

where  the  u  subscript  denotes  that  conditioning  on  the  unburned  side  of  the  flame  front  is 
being  considered.  Once  a  specific  model  for  sy  u  is  invoked,  Eq.  (6)  can  be  solved  locally  in 
an  LES  to  determine  the  value  of  a  coefficient  in  that  model  [5]. 

Both  \Vg\  and  |V{?|  are  underresolved  in  actual  LES  computations.  These  terms, 
therefore,  cannot  be  used  directly  in  Eq.  (6).  An  accurate  method  of  approximating  these 
terms  based  on  information  about  the  area  of  the  filtered  flame  front  has  been  developed 
in  the  last  two  reporting  periods  [5].  The  physics  behind  this  method  are  illustrated  in 
Fig.  1,  where  a  turbulent  flame  front  is  shown  after  having  been  filtered  with  increasingly 
wide  filter  kernels.  As  the  filter  width  increases,  the  amount  of  small  scale  structures  that 
remain  grows  smaller.  The  turbulent  burning  velocity,  therefore,  must  increase  to  ensure 
that  the  flame  front  consumes  the  correct  volume  of  unburned  gas.  This  mass  conservation 
criterion  is  accounted  for  inherently  when  filtered  flame  areas  are  used  to  represent  |V£|  in 
Eq.  (4). 
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Figure  1:  Filtered  realizations  of  a  turbulent  flame  front  propagating  at  a  speed  of  2.5  times  the  integral 
scale  velocity  fluctuation.  The  realizations  are  filtered  using  kernels  with  widths  of,  from  left  to  right,  A/t] 
=  0,  16,  32,  and  64. 


1.2.2  Turbulent  Burning  Velocity  Model  Performance 


Figure  2:  Snapshots  from  a  DNS  of  front  prop¬ 
agation.  The  level  set  is  the  wrinkled  surface, 
and  the  cut  plane  shows  voriticity  magnitude. 


A  direct  numerical  simulation  (DNS)  of  a  front  propagating  in  forced 
homogeneous  isotropic  turbulence  was  performed 
to  validate  the  dynamic  turbulent  burning  ve¬ 
locity  model.  The  DNS  was  computed  on  a 
512x256x256  mesh  at  a  Reynolds  number  of 
Re\  =  48.  The  simulation  was  run  at  a  constant 
density  to  isolate  the  propagation  model  from  gas 
expansion  effects.  A  level  set  was  used  to  describe 
the  flame  front,  and  a  laminar  burning  velocity  of 
approximately  one- third  of  the  maximum  velocity 
fluctuation  magnitude  was  prescribed.  An  image 
from  this  simulation  is  shown  in  Figure  2. 

The  DNS  was  used  to  evaluate  the  filtered 
flame  area  model  for  the  |V£|  term  in  Eq.  (4). 

The  left  hand  side  of  Fig.  3  shows  the  time  evo¬ 
lution  of  both  the  model  and  |V£|  at  a  variety  of 

filter  resolutions.  This  plot  demonstrates  that  the  model  accurately  captures  the  gradients 
of  the  Q  variable  as  a  function  of  time.  The  model  performs  well  at  all  filter  resolutions  and 
correctly  describes  how  the  gradient  fluctuations  scale  with  filter  width. 

The  right-hand  side  of  Fig.  3  shows  the  performance  of  the  proposed  dynamic  turbulent 
burning  velocity  model  at  three  canonical  filter  combinations.  The  filter  width  A/77  =  5  is 
the  case  that  is  most  representative  of  an  LES,  and  in  this  case  the  dynamic  model  performs 
much  better  than  the  model  in  which  a  constant  coefficient  is  used.  The  A/77  —  2  case  does 
not  represent  a  realistic  LES,  but  it  does  serve  to  verify  that  the  model  correctly  asymptotes 
to  a  fully  resolved  limit  when  the  filter  width  becomes  small.  Many  constant  coefficient  LES 
models  do  not  correctly  describe  the  burning  velocity  as  this  limit  is  approached.  Finally, 
the  case  in  which  A/77  —  256  demonstrates  that  the  model  performance  break  downs  for 
certain  sets  of  conditions.  These  conditions  occur  when  the  assumed  form  of  the  Sfu 
expression  that  is  used  in  Eq.  (6)  is  overtaxed.  This  limiting  case  will  need  to  be  addressed 
in  future  work. 
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Figure  3.  Left:  Time  history  of  |V£|  (—  -)  and  of  the  filtered  flame  area  per  filter  volume  Aa/Va  ( - ) 

in  the  DNS.  From  top  (black)  to  bottom  (green),  the  filter  width  is  A/tj  =  4,  8,  16,  and  32.  Right:  Time 
history  of  s^u,  from:  DNS  data  (oo),  Constant  coefficient  model  («),  Dynamic  model  with:  A/rj  =  2  and 
A'/t]  =  512  (  ••  -),  A/r]  =  5  and  A'/v  =  10  ( - ),  Afrj  =  256  and  A'/t)  =  512  ( - ). 


The  dynamic  turbulent  burning  velocity 
model  also  was  applied  in  an  LES  of  a  piloted  tur¬ 
bulent  bunsen  flame.  A  schematic  of  this  LES  is 
shown  in  Fig.  4.  The  contour  cut  plane  shows  the 
temperature  field,  while  the  isocontour  shows  the 
level  set  that  represents  the  flame  front.  The  level 
set  is  colored  by  the  ratio  of  the  filtered  and  test 
filtered  flame  front  areas  that  are  used  to  compute 
the  |V£?|  terms  in  Eqs.  (4)  and  (6).  This  ratio  is 
spatially  dependent,  and  this  dependence  is  in¬ 
troduced  to  the  flame  speed  through  the  dynamic 
coefficient  calculation. 

Figure  5  shows  that  the  use  of  the  dynamic 
procedure  strongly  affects  the  predicted  burning 
velocities.  The  left  plot  in  the  figure  shows  the 
predicted  LES  burning  velocity  when  a  constant  eas 
model  coefficient  (a)  is  used.  The  right  plot  shows 

the  predicted  burning  velocity  when  the  coefficient  is  calculated  using  a  dynamic  procedure. 
In  this  LES,  the  dynamic  model  reduces  the  magnitude  of  the  turbulent  burning  velocity 
by  as  much  as  50%  of  the  laminar  burning  velocity  sl,u.  This  reduction  can  be  attributed 
to  the  dynamic  model’s  recognition  of  the  relatively  high  resolution  of  the  LES.  The  dif¬ 
ferences  between  the  dynamic  and  constant  coefficient  models  emphasize  the  importance  of 
accounting  for  local  filtered  physics  in  LES  combustion  models. 


Figure  4:  Instantaneous  snapshot  from  an  LES 
of  the  F3  flame.  The  cut  plane  shows  a  contour 
plot  of  temperature,  which  ranges  from  300  K 
(black)  to  2200  K  (bright  yellow).  The  pre¬ 
mixed  flame  front  is  colored  according  to  the 
ratio  of  the  filtered  and  test  filtered  flame  ar- 


8 


I 

I  Figure  5:  Instantaneous  scatter  plots  of  the  turbulent  burning  velocity  as  a  function  of  turbulent  velocity 

fluctuations  in  the  bunsen  flame  LES.  Left:  is  computed  using  a  constant  value  of  the  model  coefficient 

a.  Right:  u  is  computed  locally  in  the  LES  using  a  dynamically  determined  value  of  a. 


1.2.3  Flame  Structure  Model 

Level  set  models  are  advantageous  because  they  describe  flame  front  propagation  at  arbi¬ 
trary  mesh  resolutions  without  altering  how  turbulence  affects  the  front,  but  they  have  one 
i  inherent  drawback.  Because  a  level  set  represents  a  2-D  surface,  it  cannot  be  used  to  de¬ 

scribe  the  interior  structure  of  a  premixed  flame.  A  field  quantity  such  as  a  progress  variable 
also  must  be  considered  if  this  interior  structure  is  to  be  described  accurately,  but  the  only 
way  to  protect  a  progress  variable  equation  from  numerical  errors  in  premixed  combustion 
is  to  couple  it  to  a  level  set  at  the  flame  front.  In  response  to  this  problem,  a  coupling 
procedure  was  developed  that  can  describe  3-D  flame  structure  and  account  correctly  for 
Hame  propagation  using  level  sets. 

The  level  set  and  progress  variable  coupling  is  performed  by  making  the  source  term 
in  a  progress  variable  (C)  transport  equation  dependent  on  the  location  of  a  level  set  [11]. 
This  dependency  is  introduced  by  assuming  that  the  unfiltered  C  variable  can  be  written 
as  a  presumed  function  of  a  distance  coordinate  (F  in  this  context),  and  that  Coe  can  be 
written  as  a  presumed  function  of  C.  A  1-D  flame  structure  is  assumed  additionally,  and 
the  uc{C(F))  function  then  can  be  integrated  against  a  filter  kernel  to  find  an  analytical 
expression, 

wc(C,  F,  A)  =  I"  Ltc\ F  =  F,C)P  (F  x,  L  A)  dF,  (7) 

J — oc  ' 


where  (ujc  |F  =  F,Cj  is  the  conditional  average  of  ujc  along  a  surface  that  is  a  distance 

F  =  F  away  from  the  filtered  flame  front  position.  P(F,x,£,  A)  is  the  PDF  that  describes 
the  likelihood  of  the  location  x  being  a  distance  F  away  from  the  front.  Figure  6  shows  a 
schematic  of  the  new  model. 

The  PDF  in  Eq.  (7)  couples  the  progress  variable  to  the  level  set  and  is  easy  to  compute. 
The  conditioned  source  term  in  Eq.  (7)  represents  more  of  a  challenge  but  can  be  described 
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Figure  6:  The  flame  structure  model. 


using  a  variety  of  chemical  models.  This  approach  led  to  excellent  predictions  of  the  behavior 
of  a  bluff-bodv- stabilized  premixed  flame  [11  .  The  flame  burned  in  the  thin  reaction  zones 
regime,  which  was  an  appropriate  validation  test  since,  in  this  regime,  turbulence  tends  to 
thicken  premixed  flame  structures.  In  future  work,  this  model  will  be  developed  further 
to  work  in  conjunction  with  standard  tabulated  chemistry  models  such  as  the  Flamelet 
Progress  Variable  approach. 

1.3  Partially  Premixed  Turbulent  Combustion 

Although  premixed  and  non-premixed  flamelet  solutions  can  be  mapped  cheaply  and  effec¬ 
tively  into  turbulent  flows,  these  mapping  approaches  become  problematic  when  partially 
premixed  or  multi-regime  combustion  environments  are  encountered  In  these  environments, 
it  becomes  unclear  whether  a  premixed,  a  non-premixed,  or  an  auto-ignition  chemical  so¬ 
lution  should  be  used  to  describe  chemistry.  For  example,  when  fuel  is  injected  into  a 
device  by  means  of  a  liquid  spray,  local  evaporation  and  mixing  rates  might  lead  to  either 
premixed  or  non-premixed  conditions.  Consequently,  methods  of  dealing  with  multi-regime 
and  partially  premixed  combustion  in  flamelet-based  LES  should  be  considered. 

The  flame  index  proposed  bv  Yamashita  et  al  [12]  is  one  of  the  only  available  indicators 
from  the  literature  that  can  be  used  to  determine  combustion  regimes.  The  flame  index 
simply  examines  whether  or  not  gradients  of  oxidizer  and  fuel  align  and  then  sets  the  local 
regime  accordingly  [12,  13.  14].  This  index  is  based,  however,  on  1-D  geometric  arguments, 
and  the  extent  to  which  it  can  be  applied  to  fully  3-D  turbulent  flames  is  not  yet  known. 
For  example,  in  certain  flow  settings  gradients  of  fuel  and  oxidizer  do  align,  but  combustion 
nonetheless  remains  diffusion-controlled.  The  flame  index  does  not  capture  the  physics  that 
are  relevant  in  these  regimes  [14]. 

1.3.1  Combustion  Regime  Indicator  Formulation 

A  new  Combustion  Regime  Indicator  (CRI)  model  has  been  developed  in  an  effort  to  capture 
more  of  the  physics  that  determine  local  burning  modes.  The  indicator  itself  can  be  derived 
using  a  generalized  flamelet-tvpe  transformation  of  variables.  The  mixture  fraction,  Z,  is 
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Figure  7:  Instantaneous  fields  from  the  swirl  burner.  Left:  Typical  LSB  experimental  image,  from  Cheng  [15, 
16];  Middle:  Three  C  isocontours  from  the  LES;  Right:  LES  level  set  representation  of  the  flame  front,  colored 
by  Z\ 


retained  as  one  of  the  transformation  coordinates  because  it  describes  the  asymptotic  non- 
premixed  regime.  The  key  insight  provided  by  the  indicator  lies  in  the  second  transformation 
coordinate,  which  must  describe  combustion  in  the  asymptotic  premixed  regime  and  be 
statistically  independent  of  Z .  A  flamelet  variable  A  that  satisfies  these  requirements  can 
be  defined  as  the  value  of  the  progress  variable  C  that  occurs  at  a  stoichiometric  mixture 
fraction  ( Z  =  Zst)  on  a  given  fiainelet.  When  the  scalar  transport  equation  for  a  reacting 
chemical  species  is  transformed  into  Z  and  A  space,  the  transformed  terms  can  be  grouped 
according  to  the  asymptotic  regime  in  which  they  appear.  The  transformation  takes  the 
form 


(  pdTC  +  d\C  [pdt A  +  (pu  -  puSL,un)  •  VA]) 

+  ( \C  [puSL'U  |VA|  -  V  •  (PVCV A)]  -  PYdlC)2 

+  iC,)3  =  d>c,  (8) 

where  group  1  terms  describe  unsteady  combustion,  group  2  terms  describe  premixed  com¬ 
bustion,  and  group  3  terms  describe  non-premixed  combustion.  These  terms  then  can  be 
evaluated  in  a  simulation  and  compared  to  one  another  to  determine  which  asymptotic 
group  of  terms  balances  a  local  chemical  source  term.  This  balance  indicates  the  regime  in 
which  a  fiame  burns. 

1.3.2  Model  Application:  LES  of  a  Low  Swirl  Burner 

The  CRI  has  been  tested  in  an  LES  of  the  low  swirl  burner  shown  in  Fig.  7.  This  burner 
has  been  the  subject  of  previous  experimental  and  computational  investigations  [4,  16, 
17]  and  is  particularly  interesting  to  investigate  from  the  standpoint  of  regime  prediction. 
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Regime  predictions  axe  interesting  because  the  cofiowing  air  in  the  burner  mixes  with  the 
primary  premixed  fuel  stream  and  changes  the  nature  of  the  combustion  regime  along 
the  downstream  direction.  Figure  7  shows  two  images  from  the  LES,  in  which  the  CRI 
has  been  used  actively  to  apply  asymptotic  combustion  solutions  [1].  In  particular,  the 
right  plot  indicates  how  the  coflowing  air  leans  out  the  primary  fuel  stream.  The  level  set 
isocontour  that  is  shown  in  this  plot  is  colored  by  mixture  fraction,  which  steadily  decreases 
between  the  nozzle  and  the  burner  exit.  The  left  two  plots  in  this  figure  show  that  the  flame 
predicted  by  the  CRI  methodology  instantaneously  compares  well  to  a  typical  experimental 
visualization  of  the  burner. 

Figure  8  shows  time-averaged  profiles  from  the  centerline  of  the  burner.  The  mean  axial 
velocity  is  shown  along  with  the  mean  and  rms  temperatures  as  a  function  of  the  burner 
axial  coordinate.  This  figure  is  particularly  noteworthy  because  it  demonstrates  that  the 
central  flow  feature  of  the  swirl  burner  is  captured  accurately  by  the  LES.  This  feature  is 
the  linearly  decreasing  centerline  axial  velocity  profile.  The  decrease  in  axial  velocity  is 
created  by  the  recirculation  region  that  the  swirling  flow  sets  up.  The  decreasing  velocity 
profile  acts  to  stabilize  the  position  of  the  burner’s  premixed  flame.  The  figure  shows  that 
the  LES  temperature  profiles  are  in  especially  good  agreement  with  the  experiments.  The 
agreement  suggests  that  the  combustion  model  is  performing  well,  and  that  the  CRI  and 
generalized  flamelet  transformation  theory  can  be  applied  consistently  in  LES. 


Figure  8:  Time  averaged  swirl  burner  centerline  profiles.  LES  solutions  ( - );  Experimentally  measured 

profiles  (o  o  o)  [17].  Axial  velocity,  axial  RMS  velocity,  temperature,  and  RMS  temperature  are  shown  as  a 
function  of  downstream  displacement  along  the  burner  centerline. 
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2  Development  of  an  Interactive  Platform  for  Generation, 
Comparison,  and  Evaluation  of  Kinetic  Models  for  JP-8 
Surrogate  Fuels 

2.1  Introduction 

Cost  and  efficiency  drive  the  design  of  combustion  devices  to  rely  more  and  more  on  numeri¬ 
cal  simulations.  As  the  methods  for  computational  fluid  dynamics  (CFD)  progress,  complex 
problems  such  as  the  simulation  of  chemically  reactive  flows  in  engines  become  tractable. 
Of  interest,  for  instance,  is  the  capability  to  predict  accurately  pollutant  emissions  from 
engines,  such  as  particulate  matter,  carbon  monoxide  CO,  and  oxides  of  nitrogen  NOXl  all 
contributing  at  different  levels  to  the  greenhouse  effect  and  global  warming,  smog,  ground- 
level  ozone,  stratospheric  ozone  depletion,  and  a  myriad  of  related  health  problems.  To 
predict  and  control  these  emissions,  the  understanding  and  the  accurate  modeling  of  chem¬ 
istry  is  tremendously  important.  Currently,  petroleum-derived  fuels  make  up  a  very  large 
portion  of  our  energy  resources,  and  their  combustion  proceeds  through  complex  highly 
non-linear  processes  involving  hundreds  of  different  chemical  compounds;  therefore,  a  de¬ 
tailed  chemical  modeling  of  real  hydrocarbon  fuels  is  excluded,  and  the  fuel  representation 
needs  to  be  simplified  drastically  to  be  included  in  numerical  simulations  of  combustion 
devices. 

A  first  stage  of  simplification  consists  of  approximating  the  fuel  by  a  well-defined  mixture 
of  a  few  components  (surrogate  fuels)  that  will  match  some  physical  or  chemical  properties 
of  the  real  fuel.  Using  surrogate  fuels  in  lieu  of  real  fuels  presents  numerous  advantages, 
among  which  are  the  reproducibility  of  experiments  and  the  possibility  of  formulating  chem¬ 
ical  models  suitable  for  CFD,  but  even  with  this  simplification,  deriving  chemical  models  for 
surrogate  fuels  that  can  be  used  in  CFD  simulations  remains  a  real  challenge  for  many  rea¬ 
sons.  First,  hydrocarbon  fuels  obtained  through  crude  oil  refining  processes  are  required  to 
satisfy  a  number  of  physical  or  chemical  criteria  that  sometimes  are  formulated  only  loosely; 
therefore,  their  composition  can  vary  significantly  among  feed- stocks,  and  at  best,  only  av¬ 
erage  fuel  properties  are  known.  Then  comes  the  question  of  how  to  define  the  surrogate 
compositions.  Guidelines  and  targets  have  to  be  developed  to  select  appropriate  individual 
components  and  their  respective  contribution  to  the  surrogate  mixture.  Chemical  modeling 
for  these  surrogates  is  a  major  challenge  as  well.  State-of-the-art  detailed  kinetic  models 
can  comprise  on  the  order  of  a  thousand  different  species  and  several  thousand  reactions, 
even  for  single  components.  Uncertainties  in  kinetic  data  mean  that  chemical  mechanisms 
from  different  sources  will  represent  similar  reaction  pathways  differently,  rendering  the 
merging  of  detailed  mechanisms  to  create  multi- component  mechanisms  extremely  difficult. 
Additional  complications  come  from  the  fact  that  chemical  modeling  of  single  components 
and  fuel  surrogate  compositions  likely  will  evolve,  and  previously  derived  models  will  be¬ 
come  obsolete.  Finally,  current  computational  resources  and  numerical  combustion  models 
put  severe  restriction  of  the  size  and  detail  allowable  for  the  chemistry.  Detailed  chemical 
models  must  be  reduced  considerably  in  size  before  their  implementation  in  CFD  codes  can 
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be  considered. 


2.1.1  Review  of  Objectives 

The  present  effort  aims  to  provide  some  solutions  to  the  issues  raised  above.  The  work  is 
organized  around  two  major  axes:  an  efficient  chemical  reduction  strategy  and  a  systematic, 
modular  approach  to  derive  surrogate  chemical  models.  The  former  is  crucial  to  incorpo¬ 
rate  a  certain  level  of  detail  in  the  combustion  models  of  CFD  codes,  while  flexibility  and 
consistency  appear  essential  to  go  beyond  the  current  empirical  stage  of  the  surrogate  fuel 
approach. 

A  previous  AFOSR-funded  project  led  to  the  development  of  the  Directed  Relation 
Graph  with  Error  Propagation  (DRGEP)  reduction  method,  which  automatically  eliminates 
from  a  kinetic  mechanism  the  species  and  reactions  that  do  not  contribute  significantly  to 
the  overall  dynamics  of  a  chemical  process;  however,  this  single  approach  showed  limitations, 
especially  in  the  case  of  kinetic  mechanisms  involving  a  large  number  of  branched  species. 
The  first  objective  of  the  present  work,  therefore,  is  to  improve  the  DRGEP  method  and 
complement  it  with  additional  reduction  strategies  such  as  a  chemical  lumping  method  and 
the  automatic  introduction  of  simplifying  assumptions  such  as  quasi-steady  states. 

The  next  objective  is  the  development  of  a  modular  framework,  called  the  component 
library  approach,  to  formulate  surrogate  compositions  optimized  for  a  given  application  and 
specific  hydrocarbon  fuel  automatically,  and  generate  validated  multi-component  reduced 
models  for  these  surrogates.  Emphasis  is  placed  on  developing  reduced  models  for  JP-8. 
To  achieve  this  goal,  a  preliminary  version  of  the  component  library  developed  under  our 
previous  AFOSR  grant  needed  to  be  extended  to  include  compounds  relevant  to  jet  fuels, 
such  as  long-chain  alkanes,  cyclo- alkanes,  and  aromatic  components.  To  obtain  the  small¬ 
est  possible  kinetic  models,  the  potential  of  a  modular  approach  was  explored,  in  which 
modules  describing  for  example  low- temperature  ignition  or  soot  formation  were  developed 
and  included  in  the  surrogate  kinetic  models  only  if  necessary.  Finally,  the  viability  of  the 
overall  surrogate  modeling  strategy  needed  to  be  assessed  through  the  development  of  a 
reduced  jet  fuel  surrogate  model. 


2.2  Automatic  Reduction  of  Detailed  Chemical  Kinetic  Mechanisms 

Three  different  reduction  techniques  have  been  developed,  each  of  them  addressing  a  dif¬ 
ferent  aspect  of  kinetic  reduction.  The  first  one,  called  Directed  Relation  Graph  with 
Error  Propagation  (DRGEP)  method,  eliminates  species  and  reactions  [18],  the  second  one 
compacts  information  through  chemical  isomer  lumping  [19],  and  the  third  one  introduces 
quasi-steady  state  assumptions  using  a  criterion  based  on  lifetime  analysis  and  DRGEP  [20]. 
The  different  techniques  are  combined  into  an  integrated  approach  that  is  shown  to  reduce 
drastically  the  size  of  large-scale  kinetic  mechanisms. 
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2.2.1  Elimination  Stage:  The  Directed  Relation  Graph  with  Error  Propagation 
Method 

The  goal  of  the  elimination  stage  is  to  identify,  for  any  number  of  species  in  the  skeletal 
mechanism,  Arskeh  a  group  of  species  of  size  N ^  =  A^et  -  A^l  that  can  be  removed  with 
minimal  impact  on  the  targets.  The  choice  of  these  species  is  done  here  by  defining  appro¬ 
priate  importance  coefficients  for  each  species  based  on  the  production  and  consumption 
rates,  which  are  evaluated  using  results  obtained  from  the  detailed  mechanism.  The  species 
with  the  ATrm  lowest  importance  coefficients  are  removed  from  the  mechanism,  and  a  skeletal 
mechanism  of  size  Ns^e\  is  created  by  removing  from  the  detailed  mechanism  any  reaction 
in  which  a  removed  species  appears  as  a  reactant  or  as  a  product. 

Direct  interaction  coefficients  Direct  interaction  coefficients  are  defined  as  the  measure 
of  the  coupling  between  two  species  that  are  related  directly  through  an  elementary  reaction, 
that  is,  two  species  that  appear  concurrently  in  the  same  reaction.  In  the  DRGEP  method, 
the  coupling  coefficient  between  two  directly  related  species  A  and  B  is  estimated  as  follows: 


max  (Pa,  Ca) 


(9) 


where 

Pa  =  ^  max  (0,  ViyAVt),  and  Ca  =  ^  max  (0,  —  (10) 

i=l,  rift  i=l,rifl 

Equation  (9)  provides  an  estimate  of  the  impact  that  removing  one  species  has  on 
the  calculated  concentration  of  the  remaining  species;  however,  the  goal  of  the  reduction 
procedure  is  to  remove  the  largest  possible  set  of  species  from  the  mechanism  while  keeping 
errors  below  a  given  tolerance.  Considering  one  species  independently  of  the  group  of 
removed  species  in  which  it  will  belong  eventually  might  lead  to  a  very  inaccurate  estimate 
of  the  importance  of  each  species.  This  observation  leads  to  the  extension  of  Eq.  (9)  given 
a  set  of  removed  species: 


rAB,{S} 


Sz=l,nH  ui.A^i^{S} 

ma x(PA,CA) 


(11) 


where  {<S}  is  the  set  of  species  already  removed.  5lB  is  unity,  if  the  ith  reaction  involves 
B  or  any  species  in  subset  {<S},  and  0  otherwise. 


Error  Propagation  The  goal  of  the  elimination  stage  is  to  remove  species  and  reactions 
that  negligibly  impact  a  user-defined  set  of  targets.  Intuitively,  the  farther  away  from  the 
target  a  species  is,  the  smaller  the  effect  of  changing  or  removing  this  species  should  be. 
To  reflect  this  error  propagation  effect,  a  geometric  damping  has  been  introduced  in  the 
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selection  procedure;  therefore,  we  define  a  path-dependent  coefficient: 


n—  1 

rAB,P=UrSiSi+,  (12) 

i=  1 

and  a  global  importance  coefficient 

Rab=  max  rAB,p.  (13) 

all  paths  p 

If  some  error  is  introduced  in  the  prediction  of  a  species  B ,  the  longer  the  way  this  error  has 
to  propagate  to  reach  the  target  A ,  the  smaller  its  effect  will  be  typically.  This  technique  is 
target-oriented  and  is  expected  to  provide  a  finer  selection  of  the  chemical  paths  necessary 
for  the  accurate  prediction  of  the  set  of  targets  by  keeping  species  associated  with  large  R 
coefficients  and  discarding  species  with  small  R  coefficients.  A  similar  procedure  is  adopted 
to  reduce  the  number  of  reactions  [18]. 

Integrity  Check  Every  intermediate  species  in  a  skeletal  mechanism  must  have  at  least 
one  production  and  one  consumption  path.  During  the  DRGEP  reduction  process,  some 
species  might  fail  this  requirement,  especially  for  high  reduction  ratios  in  the  context  of 
complex,  highly  non-linear  kinetic  schemes.  Two  basic  observations  can  be  made.  In  a 
closed  system,  any  intermediate  species  that  is  not  produced  anymore  remains  at  its  initial 
zero  concentration  and  can  be  removed  from  the  mechanism.  On  the  other  hand,  if  a 
species  is  produced,  but  not  consumed  anymore,  it  creates  a  sink  of  mass  that  may  impact 
greatly  the  final  concentration  of  the  products.  This  complex  non-linear  behavior  cannot 
be  detected  by  a  method  based  solely  on  the  analysis  of  the  detailed  production  rates; 
therefore,  a  simple  algorithm  has  been  designed  to  prevent  these  situations  from  occurring. 
A  list  of  species,  sorted  by  order  of  importance  for  the  targets,  is  obt  ained  first  by  computing 
the  DRGEP  coefficients  using  Eq.  (9).  Then  this  list  is  modified  slightly  so  that,  for  any 
value  of  the  cutoff  parameter,  the  group  of  species  kept  in  the  skeletal  mechanism  forms  a 
consistent  chemical  scheme  with  no  truncated  paths. 

Validation  Before  applying  the  above  reduction  methodology,  the  validity  of  the  error 
propagation  assumption  in  the  DRGEP  method  needs  to  be  appraised.  This  assumption 
states  that  the  effect  on  a  target  introduced  by  the  removal  of  a  species  can  be  approxi¬ 
mated  through  geometric  damping  along  the  directed  relation  graph,  from  the  target  to  the 
removed  species,  which  can  be  written  as 

EtocRtSi  (14) 

where  Et  is  the  error  between  the  prediction  of  the  target  using  the  detailed  and  skeletal 
mechanisms  and  Rts  measures  the  importance  of  a  species  S  with  respect  to  the  target,  as 
defined  by  Eqs.  (9)  and  (13).  This  proportionality  can  be  verified  a  posteriori  using  a  prac¬ 
tical  case.  The  coefficients  Rts  were  computed  for  the  adiabatic,  isochoric  auto-ignition  of 
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a  stoichiometric  mixture  of  iso-octane  and  air  at  13  bar  and  1000  K,  with  iso-octane  as  the 
only  target.  The  mechanism  used  for  the  simulation  was  the  iso-octane  oxidation  scheme 
of  Curran  et  al.  [21].  Figure  9  shows  the  correlation  between  this  error  and  the  computed 
coefficients  obtained  using  error  propagation.  For  comparison,  the  correlation  between  the 
error  and  the  coefficients  obtained  without  error  propagation  is  shown  also  in  the  same 
figure.  The  solid  line  represents  the  optimal  case,  that  is,  a  hypothetical  parameter  whose 
value  would  be  exactly  equal  to  the  error  introduced  if  the  species  were  removed.  The 


Figure  9:  Correlation  between  the  error  introduced  in  the  fuel  prediction  and  the  species  coefficients  Ricg  s 
obtained  with  and  without  error  propagation  during  the  isochor,  adiabatic  auto-ignition  of  a  stoichiometric 
mixture  of  iso-octane  and  air  at  13  bar  and  1000  K. 

errors  introduced  by  removing  individual  species  correlate  extremely  well  with  the  error 
propagation  coefficients,  with  a  small  scatter  in  the  data,  whereas  the  correlation  is  not  as 
obvious  without  error  propagation.  Figure  9  also  shows  that  the  error  propagation  method 
leads  to  an  order  unity  coefficient  in  Eq.  (14).  This  observation  means  that  the  impor¬ 
tance  coefficients  evaluated  by  DRGEP  are  a  direct  measure  of  the  error  in  the  resulting 
mechanism. 

The  efficiency  of  the  group-based  coefficients  and  the  integrity  check  are  illustrated 
next.  The  detailed  mechanism  for  iso-octane  oxidation  from  Curran  et  al.  [21]  is  reduced 
for  a  single  initial  physical  condition,  the  homogeneous,  adiabatic  auto-ignition  at  constant 
volume  of  a  stoichiometric  mixture  of  iso-octane  and  air  at  an  initial  pressure  of  13  bar  and 
an  initial  temperature  of  625  K.  The  targets  are  fuel,  CO,  COo,  and  temperature.  Figure  10 
shows  the  evolution  of  the  error  in  ignition  delay  time  and  in  the  final  mass  fraction  of  carbon 
monoxide  as  functions  of  the  number  of  species  kept  in  the  skeletal  mechanism.  When  the 
definition  of  Eq.  (9)  is  used,  the  error  in  the  final  mass  fraction  of  CO  in  the  system  quickly 
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Figure  10:  lso-octane  auto-ignition  at  low  temperature.  Evolution  of  the  error  in  ignition  delay  time  (lines 
with  symbols)  and  final  mass  fraction  of  CO  (plain  lines)  when  the  group-based  coefficients  and  the  integrity 
check  algorithm  are  used  (open  symbols,  dashed  lines)  or  not  (filled  symbols,  solid  lines). 


reaches  a  few  percent  and  keeps  increasing.  This  behavior  is  due  to  truncated  chemical 
paths  appearing  as  species  axe  removed.  Carbon  mass  accumulates  in  large  quantities 
in  intermediate  species,  which  shifts  considerably  the  chemical  equilibrium  of  the  system. 
When  group-based  coefficients  and  integrity  check  are  included  in  the  reduction  process 
(the  coefficients  are  recomputed  once  every  50  species  removed),  the  error  in  the  final  mass 
fraction  of  CO  remains  extremely  small,  and  the  error  in  ignition  delay  time  is  improved 
considerably. 

To  demonstrate  the  full  capabilities  of  the  DRGEP  method,  the  same  very  large  mecha¬ 
nism  for  isooctane  oxidation  [21]  is  reduced  for  adiabatic  autoignition  at  constant  volume 
over  a  large  range  of  initial  conditions  relevant  for  engine-related  applications  (ignition  de¬ 
lay  times  less  than  one  second).  The  initial  conditions  include  equivalence  ratios  between 
0.5  and  2,  pressures  between  1  bar  and  40 bar,  and  temperatures  between  600  K  and  1500  K. 
The  detailed  mechanism  comprises  850  species  and  7212  reactions.  Targets  for  the  reduction 
are  fuel  ir ChIIis,  major  products  CO  and  CO2,  and  temperature.  The  DRGEP  coefficients 
for  temperature  are  evaluated  using  heat  release  data.  Overall,  errors  are  shown  to  in¬ 
crease  monotonically  as  the  number  of  species  is  reduced.  The  total  number  of  species  kept 
in  the  skeletal  mechanism  is  chosen  so  that  the  maximum  error  over  all  targets  is  about 
15%,  which  corresponds  to  196  remaining  species  and  1762  remaining  reactions,  forward 
and  backward  counted  separately.  Following  this  first  step  of  reduction  and  with  the  same 
accuracy  requirement,  additional  non-necessary  reactions  are  removed,  and  the  resulting 
skeletal  mechanism  comprises  195  species  and  802  reactions,  that  is,  a  reduction  by  a  factor 
of  4.35  of  the  number  of  species  and  a  reduction  by  a  factor  of  9  of  the  number  of  reactions. 
Table  1  provides  an  estimate  of  the  error  introduced  in  ignition  delay  times  over  the  entire 
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reduction  domain  by  this  first  stage  of  reduction. 


A  Species 

A  Reactions 

Maximum  Error  [% 

Average  Error  [%] 

196 

1762 

15.89 

6.02 

195 

802 

14.96 

5.55 

Table  1:  Maximum  and  average  errors  introduced  by  the  elimination  stage  of  reduction,  for  initial  conditions 
with  pressures  between  1  and  40  bar,  equivalence  ratios  between  0.5  and  2,  and  temperatures  between  600  K 
and  1500  K. 


2.2.2  Compacting  Stage:  A  Chemical  Lumping  Approach 

A  second  important  aspect  of  kinetic  reduction  aims  at  describing  the  system  in  terms 
of  a  reduced  number  of  variables,  called  lumped  variables,  through  a  linear  or  non-linear 
transformation.  In  the  present  work,  a  general  automatic  lumping  approach  was  derived 
that  directly  uses  simulation  results  obtained  from  the  detailed  mechanism  to  generate  a 
lumped  scheme  valid  over  a  user-specified  range  of  conditions.  The  method,  applied  here  to 
isomer  lumping  in  hydrocarbon  oxidation  kinetic  schemes,  does  not  rely  on  equilibrium  or 
quasi-  steady  state  assumptions.  In  addition,  the  resulting  lumped  mechanisms  are  suitable 
for  direct  use  in  standard  chemistry  software.  In  the  following,  the  proposed  lumping 
procedure  will  be  detailed  briefly,  and  the  quality  of  the  resulting  lumped  schemes  will  be 
assessed. 


Proposed  Modeling  Approach  Suppose  the  original  set  of  species  S  =  {Sj,j=i..jv5 }  is 
sub-divided  into  Ns  lumped  groups  C *  ^  ^  ~  .  For  each  lumped  group  ,  a  representative 
species  is  defined  as  a  linear  combination  of  the  species  of  this  group,  such  that: 


£  M  ■ 

iecf 


(15) 


We  define  the  relative  contribution  of  species  Si  to  its  group  as: 


Oi  = 


(16) 


The  reaction  rate  of  a  reaction  j  can  be  expressed  using  the  lumped  variables: 


Ns  ,  NS 
i=l  1=1 


1*3 


(17) 


where  At  this  stage,  the  reactions  that  have  become  identical  after  the 

lumping  of  species  also  can  be  lumped  together,  the  rate  of  each  lumped  reaction  being 
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defined  as: 


11 

M 

£ 

II 

to) 

fc»i 

with  kj  =  Y  (^117^1  ' 

(18) 

/=1 

jec*}  V  *=1  / 

This  transformation  corresponds  to  an  exact  linear  lumping  if  kj  jy  ,  and,  therefore, 
the  relative  contributions  are  known  functions  of  time  and  space.  Unfortunately, 

these  functions  usually  are  unknown,  and  a  closed  form  has  to  be  assumed.  In  this  work, 
rather  than  specifying  the  relative  distributions  as  functions  of  physical  space  and  time, 
models  for  these  distributions  will  be  formulated  in  terms  of  the  state  space: 


T,P, 


5/ 


1=1. .Ns 


(19) 


Suppose  that  a  subset  II  of  the  state  space  variables  is  chosen  to  parameterize  a  set 
of  data  a.  In  our  case,  rr  refers  to  the  relative  distributions  OLij=\..Ns  or  to  the  lumped 
rates  ^  J=1  ^  In  most  cases,  no  explicit  relationship  between  the  data  and  the  chosen 
state  parameters  can  be  found;  therefore,  a  model  fn  (Ft)  has  to  be  formulated  that  will 
approximate  the  actual  data  by  a  function  of  the  parameters  included  in  II.  Both  the  choice 
of  the  set  of  parameters  and  the  definition  of  the  function  /  will  impact  the  quality  of  the 
model. 

A  promising  approach  is  to  use  non-linear  data  modeling  tools  such  as  artificial  neural 
networks,  which  can  be  seen  as  parametric  functions  whose  weights  are  adjusted  by  training 
the  network  to  minimize  errors  (22).  This  technique  retains  the  dependence  of  the  data  on  a 
large  number  of  state  variables,  thus  introducing  very  little  error  in  the  model;  however,  the 
training  process  is  not  an  easy  and  straightforward  procedure,  and  the  resulting  mechanism 
would  not  be  used  easily  in  standard  chemistry  packages,  such  as  the  Chemkin  libraries, 
as  extra  routines  need  to  be  provided  to  evaluate  the  lumped  rate  constants;  therefore, 
for  practical  reasons,  the  approach  will  be  demonstrated  and  analyzed  by  considering  only 
constant  values  or  temperature  dependent  values  of  az.  To  increase  accuracy,  the  modeling 
procedure  will  be  applied  directly  to  the  lumped  rate  coefficients  kj  instead  of  the  relative 
contributions  a*.  Then,  the  lumped  rate  coefficients  can  be  expressed  in  Arrhenius  form 
and  the  lumped  reactions  incorporated  readily  into  the  lumped  mechanism  that  retains  the 
same  format  as  the  detailed  mechanism.  Accordingly,  kj  will  be  modeled  using  a  basis 
function  of  the  form: 

fkj  (T)  =  pT^e-zr  .  (20) 

The  determination  of  the  unknown  coefficients  /?,  7,  and  5  is  done  using  the  detailed  sim¬ 
ulation  data.  To  illustrate  the  soundness  of  this  simplifying  assumption,  the  data  for  the 
iso-octyl  isomers  obtained  using  the  mechanism  of  Curran  et  al.  [21]  are  projected  onto  the 
temperature  and  shown  in  Fig.  11.  The  role  played  by  these  isomers  in  the  overall  dynamics 
of  the  system  might  be  more  or  less  important  depending  on  the  configurations.  A  global 
interaction  coefficient  was  defined  as  part  of  the  DRGEP  methodology  presented  in  [23] 
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that  quantifies  the  impact  of  a  species  on  major  parameters  such  as  ignition  delay  time, 
laminar  burning  velocities,  or  fuel  consumption.  This  DRGEP  coefficient  R  is  computed  for 
all  cases  and  used  as  a  weighting  factor  to  analyze  the  relative  contribution  of  each  iso-octyl 
isomer.  In  Fig.  11,  darker  colors  mean  larger  coefficients  /?,  that  is,  a  higher  importance  of 
the  radicals  in  the  mechanism.  These  data  demonstrate  that  keeping  only  the  temperature 
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Figure  11:  Iso-octyl  radical  distribution  as  function  of  the  temperature  for  homogeneous  and  flame  con¬ 
figurations.  Darker  colors  mean  higher  DRGEP  coefficients.  Dashed  lines  correspond  to  thermodynamic 
equilibrium  ratios. 

dependence  of  the  distribution  functions  a  is  a  good  approximation  at  low  temperatures, 
whereas  some  scatter  is  observed  at  higher  temperatures;  however,  less  accuracy  in  the 
isomers  distribution  will  have  little  impact  on  the  prediction  of  the  targets,  as  indicated 
by  smaller  DRGEP  coefficients.  Also  shown  in  Fig.  11  are  the  relative  contributions  com¬ 
puted  using  the  equilibrium  constants  of  the  isomerization  reactions.  Although  the  trends 
are  similar,  using  the  pseudo-equilibrium  assumption  to  evaluate  a  might  introduce  large 
inaccuracies  in  the  reaction  rates. 

To  maintain  high  accuracy  in  hydrocarbon  oxidation  mechanisms,  the  overall  organi¬ 
zation  of  the  detailed  mechanism  should  be  kept  in  the  lumped  mechanism.  For  example, 
representative  elementary  steps  of  each  specific  class  of  reaction  should  be  present  regard¬ 
less  of  the  extent  of  reduction.  The  presence  of  these  elementary  steps  allows  the  reduced 
mechanism  to  be  used  to  get  insight  into  the  dynamic  chemical  processes  occurring  during 
combustion.  A  simple  way  to  achieve  this  result,  which  is  followed  here,  is  to  lump  together 
chemical  isomers,  as  these  species  have  potentially  similar  formation  and  decomposition 
pathways.  Additionally,  for  some  groups  of  isomers,  taking  into  account  the  size  of  the  ring 
in  the  transition  states,  as  suggested  by  Ahmed  et  al.  [24],  improved  results  dramatically, 
with  a  small  loss  in  the  degree  of  reduction. 

Results  and  Validation  The  skeletal  mechanism  for  iso-octane  oxidation  obtained  through 
species  and  reaction  elimination  [23]  is  used  as  a  starting  mechanism  to  demonstrate  the 
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efficiency  of  the  lumping  approach.  This  mechanism  contains  195  species  and  802  reactions, 
backward  and  forward  reactions  being  counted  separately.  Among  those  species,  27  lumped 
groups  involving  88  chemical  isomer  species  can  be  identified,  leading  to  a  reduction  of  the 
number  of  species  from  195  to  135  if  all  groups  of  isomers  are  lumped  and  a  reduction  of 
the  number  of  reactions  to  611. 

In  addition  to  the  auto-ignition  configurations  used  to  generate  Fig.  11,  the  lumped 
mechanism  obtained  above  with  the  proposed  method  was  used  to  simulate  atmospheric  plug 
flow  reactors  modeled  by  isobaric  homogeneous  systems,  and  atmospheric  omvdimensional 
laminar  premixed  flames.  Table  2  provides  the  errors  obtained  over  all  conditions  between 
skeletal  and  lumped  schemes  for  ignition  delay  times  and  laminar  burning  velocities,  which 
are  two  parameters  describing  the  global  behavior  of  the  considered  systems.  Virtually  no 
error  on  the  final  concentration  of  the  products  was  observed.  The  small  errors  show  that  the 


Error  in  ng  [%] 


Error  in  Sl 


Max 

16.19 

Avg 

4.75 

Max 

1.05 

Avg 

0.67 

Table  2:  Comparison  of  global  parameters  between  skeletal  and  lumped  iso-octane  mechanisms. 


prediction  of  global  parameters  is  quite  accurate.  More  remarkably,  even  though  the  lumped 
rate  coefficients  were  modeled  using  homogeneous  reactor  data  only,  the  resulting  scheme 
is  applicable  to  different  configurations  such  as  flames;  however,  additional  comparisons  are 
required  to  ensure  that  the  lumping  does  not  affect  fundamental  aspects  of  the  original 
mechanism,  for  example,  mass  fluxes  between  species.  The  adequate  representation  of  a 
group  of  isomers  can  be  assessed  by  comparing  the  sum  of  the  concentration  of  isomers  of 
the  lumped  group  to  the  actual  concentration  of  the  representative  species.  For  an  accurate 
lumping,  these  two  quantities  should  be  equal,  as  indicated  by  the  fundamental  definition 
of  the  representative  species  in  Eq.  (15).  This  property  is  demonstrated  in  Fig.  12  for  the 
iso-octane  lumped  mechanism  for  an  auto-ignition  case,  which  shows  that  the  concentration 
of  the  representative  octyl  radical  compares  very  well  to  the  sum  of  the  four  octyl  radicals 
lumped  together. 

2.2.3  Introduction  of  Modeling  Assumptions:  The  Quasi-Steady  State  Ap¬ 
proximation 

To  increase  the  speed  up  of  the  skeletal  mechanism  obtained  through  species  and  reaction 
elimination  and  chemical  lumping  even  further,  a  straightforward  strategy  is  the  introduc¬ 
tion  of  quasi-steady  state  assumptions  that  replace  part  of  the  differential  equations  by 
algebraic  equations,  which  are  much  faster  to  evaluate.  In  the  present  work,  suitable  QSS 
species  are  identified  systematically  through  the  evaluation  of  a  steady-state  parameter 
based  on  lifetime  analysis  and  DRGEP  coefficients  [18].  The  steady  state  parameter  for  a 
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Figure  12:  Comparison  of  iso-octyl  isomers  and  lumped  reprasentative  species  concentrations  during  auto¬ 
ignition. 


species  S  can  be  expressed  as: 


Qs  (f)  =  ar/?T5[5]r,s' , 


(21) 


where  [S]  is  the  concentration  of  species  S  and  the  lifetime  of  this  species  can  be  expressed 
in  terms  of  its  production  and  consumption  rates 


TS  =  ~ 


\d(Ps-Cs)  1 

.  5  [5]  . 


(22) 


The  species  with  small  Q  values  for  all  cases  can  be  set  in  steady  state.  For  simplicity,  only 
linear  coupling  between  steady  state  species  is  allowed,  so  that  explicit  expressions  can  be 
written  automatically  for  direct  use  in  a  combustion  code. 

To  assess  the  validity  of  the  method,  the  values  for  the  steady  state  parameter  Q  given 
by  Eq.  (21)  were  computed  using  the  195  species  isooctane  mechanism  obtained  after 
the  elimination  stage,  for  all  cases  and  targets  used  in  the  validation  of  this  mechanism. 
Additionally,  the  actual  error  in  ignition  delay  time  obtained  when  setting  species  in  steady 
state  one  after  the  other  was  computed  for  each  configuration.  This  error  is  correlated  with 
the  corresponding  value  of  Q  for  the  species.  Results  are  presented  in  Fig.  13.  A  very 
small  error  of  the  order  of  10-5  is  introduced  systematically  that  is  due  to  numerics  and 
grid  resolution.  A  clear  trend  is  observed,  with  species  having  small  Q  values  introducing 
comparatively  smaller  error  than  species  with  large  Q  values.  For  example,  a  cut-off  value 
of  Q  =  io-12-25  identifies  correctly  80  out  of  83,  that  is,  96%  of  the  species  introducing  less 
than  0.2%  error  on  ignition  delay  time  when  set  in  steady  state. 

The  evaluation  of  the  coupled  set  of  algebraic  equations  corresponds  to  the  inversion 
of  a  matrix-vector  system.  To  optimize  the  computational  time  required  to  perform  this 
inversion,  a  search  algorithm  first  identifies  groups  of  inter-dependent  species  and  re-orders 
the  QSS  species  such  that  the  resulting  matrix  is  a  block- triangular,  easily  invertible,  matrix. 
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Figure  13:  Maximum  error  in  ignition  delay  time  against  the  steady  state  parameter  Q  of  each  species  when 
set  in  steady  state  over  the  specified  range  of  initial  conditions  used  for  the  reduction  of  the  iso-octane 
mechanism. 

2.2.4  Integration  into  a  Multi-Stage  Reduction  Strategy 

To  exploit  the  reduction  potential  of  the  techniques  presented  above  fully,  the  three  reduc¬ 
tion  methods  are  combined  through  a  multi-stage  approach  and  applied  to  reduce  the  LLNL 
comprehensive  iso-octane  oxidation  mechanism  of  Curran  et  al.  [21].  The  sample  state  space 
includes  auto-ignition,  plug-flow  reactors,  and  premixed  flames.  The  order  of  the  reduction 
techniques  is  as  follows.  First,  DRGEP  removes  negligible  species  and  reactions  until  some 
error  tolerance  of  about  15%  for  the  chosen  targets  is  reached.  Removing  more  species  at 
this  stage  would  bring  the  error  up  considerably,  as  parts  of  some  important  parallel  path¬ 
ways  would  be  removed.  Although  each  channel  individually  does  not  contribute  much  to 
the  global  fluxes,  discarding  too  many  of  them  ends  up  introducing  significant  error.  These 
pathways  are  lumped  together  in  the  next  stage  of  reduction  using  the  isomer  lumping 
method  presented  above.  Then,  an  additional  stage  of  DRGEP  can  be  applied  to  reduce 
the  number  of  species  further.  Finally,  QSSA  are  introduced  to  speed  up  the  computational 
process.  A  final  reduced  mechanism  is  obtained,  which  consists  of  57  species,  52  steady- 
state  species,  and  504  reactions.  Comparisons  of  the  results  at  various  stages  of  reduction 
are  shown  in  Figs.  14.  Overall,  the  reduced  mechanism  reproduces  correctly  the  detailed 
results  for  ignition  timing,  very  lean  oxidation  in  PFR,  and  flame  propagation.  The  errors 
introduced  by  the  reduction  are  small  considering  the  high  reduction  ratio  and  everywhere 
negligible  compared  to  the  discrepancies  with  experimental  data. 

To  appraise  the  benefits  of  the  mechanism  reduction,  average  computational  costs  are 
recorded  for  each  skeletal  or  reduced  mechanism  and  compared  to  the  cost  of  using  the 
full  detailed  model.  Results  are  shown  in  Table  3  for  homogeneous  cases.  The  observed 
computational  gain  scales  roughly  like  the  square  of  the  number  of  variables  kept  in  the 
system,  with  a  small  additional  benefit  linked  to  the  elimination  of  non-important  reac¬ 
tions.  A  detailed  flux  analysis  on  detailed  and  reduced  schemes  shows  that  the  various 
fuel  consumption  paths  and  their  relative  importance  are  retained  for  both  low  and  high 
temperature  oxidation  [20]. 
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Figure  14:  Ignition  delay  times  for  iso-octane.  Comparison  between  experiments  [25,  26,  27)  (filled  symbols), 
detailed  (solid  lines),  195  species  skeletal  (dashed  lines),  and  57  species  reduced  (open  circles)  mechanisms. 
Maximum  error  in  ignition  delay  time  is  34.94  %,  average  error  is  10.86  %.  Errors  in  the  final  concentration 
of  the  major  products  H2O  and  CO2  do  not  exceed  0.5  %. 


Number  of  species  kept 

Relative  cost 

Full  mechanism 

850 

100  % 

100 

DRGEP  (species  elimination) 

196 

23.06  % 

4.37 

DRGEP  (reaction  elimination) 

195 

22.94  % 

2.82 

Lumping 

134 

15.76  % 

2.16 

DRGEP  (extra) 

109 

12.82  % 

1.65 

QSS 

57 

6.70  % 

0.98 

Table  3:  Relative  computational  cost  for  homogeneous  simulations  of  the  intermediate  skeletal  and  final 
reduced  mechanisms  produced  during  the  integrated  multi-stage  reduction  of  a  detailed  mechanism  for  iso¬ 
octane. 
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2.3  A  Modular  Approach  to  Model  Transportation  Fuel  Chemistry 

2.3.1  Selection  of  Surrogate  Components  and  Composition 

When  designing  a  surrogate  fuel,  the  goal  is  to  identify  a  limited  number  of  hydrocarbon 
molecules  that  can  be  blended  into  useful  experimental  fuels  and  modeled  computation¬ 
ally  [28].  Numerical  considerations  imply  that  kinetic,  thermodynamic,  and  physical  data 
of  adequate  quality  axe  available  for  each  of  the  molecules  included  in  the  set.  Then,  careful 
analysis  of  the  type  of  application  for  which  the  surrogate  is  needed  must  be  done,  as  the 
application  will  determine  the  set  of  relevant  targets  the  surrogate  will  have  to  match.  Po¬ 
tential  targets  are  the  heating  value  or  energy  content  of  the  fuel,  physical  properties  of  the 
liquid  fuel,  boiling  characteristics,  reactivity,  pollutant  emissions,  and  class  composition.  A 
discussion  of  the  applicability  of  each  of  these  targets  and  how  they  can  be  evaluated  for  a 
fuel  mixture  can  be  found  in  [20]. 

Appropriate  hydrocarbon  candidates  to  be  included  in  surrogate  fuels  should  have  been 
studied  both  experimentally  and  kinetically,  so  that  validated  models  of  sufficient  accuracy 
are  available.  A  review  of  available  detailed  chemical  kinetic  models  for  the  oxidation  of 
hydrocarbon  molecules  has  been  done  by  Simmie  [29].  A  certain  number  of  such  molecules 
and  their  potential  relevance  to  the  three  major  transportation  fuels  have  been  identified  [28, 
30,  31].  Following  these  guidelines,  we  have  used  in  this  work  ?>heptane,  decane,  and  n- 
dodecane  for  the  linear  paraffin  group,  iso-octane  and  methyl-cyclohexane  as  branched 
paraffins,  and  benzene  and  toluene  as  aromatics. 

A  constrained  optimization  algorithm  has  been  developed  to  facilitate  the  formulation 
of  surrogate  composition.  The  procedure  relies  on  the  fact  that  most  properties  described 
above  are  bulk  properties  and  thus  easily  expressed  as  combinations  of  the  properties  of  the 
individual  components.  For  properties  such  as  threshold  sooting  indices,  for  which  data  may 
not  be  available  for  all  components,  a  group  contribution  approach  has  been  adopted  that 
estimates  the  properties  of  the  mixture  as  function  of  the  structural  groups  of  the  molecules 
included  in  the  mixture.  The  procedure  follows  that  described  by  Yan  et  al.  [32],  based  on 
the  initial  work  of  Benson  et  al.  [33].  Given  a  set  of  potential  components,  the  algorithm 
returns  the  optimal  composition  that  matches  best  the  user-specified  set  of  constraints. 

2.3.2  Component  Library  Approach 

Developing  a  chemical  mechanisms  for  a  transportation  fuel  surrogate  will  not  be  done 
only  once  but  will  be,  rather,  a  dynamic  process.  CFD  simulations  require  mechanisms 
that  are  as  small  as  possible.  This  requirement  definitely  sets  aside  the  uone  surrogate 
fits  all”  approach  and  favors  surrogates  tightly  tailored  to  the  considered  application,  but 
even  then,  the  composition  provided  by  the  optimization  approach  that  relies  on  global 
characteristics  may  not  be  fully  adequate  for  the  problem  at  hand,  and  further  experimental 
data  may  suggest  some  different  composition  or  additional  components.  For  example,  the 
surrogate  proposed  by  the  jet  fuel  surrogate  working  group  [31],  composed  of  50%  Tvdecane, 
25%  n-butylbenzene,  and  25%  n-butylcyclohexane,  was  tested  in  a  pressurized  flow  reactor 
and  in  a  single  cylinder  research  engine,  and  its  characteristics  in  terms  of  flame  ignition 
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and  extinction  in  counterflow  flames  were  compared  to  an  average  jet  fuel.  Although  the 
surrogate  matched  the  hydrogen-to-carbon  ratio,  it  was  consistently  much  more  reactive 
than  the  fully-blended  fuel,  and  therefore,  was  not  an  acceptable  surrogate. 

On  the  other  hand,  detailed  chemical  mechanisms  undergo  constant  modifications,  as 
the  knowledge  and  understanding  of  the  underlying  kinetic  processes  increase,  and  these 
modifications  must  be  incorporated  into  the  reduced  schemes  for  surrogate  fuels  as  well 
An  additional  challenge  lies  in  the  sizes  of  detailed  mechanisms.  If  the  detailed  kinetic 
description  of  a  single  component  involves  up  to  a  thousand  species  and  several  thousands 
reactions,  combining  several  of  these  to  form  the  model  for  a  surrogate  is  clearly  not  a 
trivial  task,  and  the  chance  of  introducing  incoherencies,  such  as  truncated  paths  or  invol¬ 
untarily  duplicated  reaction  pathways,  is  very  high.  Moreover,  mechanisms  of  this  size  are 
impractical  to  use  even  in  the  simplest  configurations.  Numerical  problems  may  arise,  and, 
most  importantly,  no  accurate  and  detailed  analysis  and  validation  of  the  mechanism  can 
be  carried  out.  Any  reduction  algorithm  is  then  rendered  inefficient  or,  at  best,  extremely 
slow. 

We  propose  here  a  methodology  to  construct  reduced  kinetic  models  for  surrogate  fuels 
that  is  based  on  two  principles: 

•  Simplification:  Each  type  of  mechanism  manipulation  has  to  be  done  at  its  simplest 
level.  For  example,  reduction  should  be  done  on  single  component  mechanisms,  as  it 
reduces  both  the  size  of  the  initial  scheme  and  the  extent  of  the  validation  domain. 
Also,  combining  kinetic  data  of  single  components  to  form  mechanisms  that  can  han¬ 
dle  mixtures  should  be  done  with  the  smallest  possible  starting  mechanisms,  ideally 
skeletal  ones. 

•  Flexibility:  The  addition  of  another  component  in  the  definition  of  the  surrogate 
mixture  should  not  require  the  entire  process  to  be  repeated.  Instead,  only  the  parts 
relevant  to  this  extra  molecule  should  be  added.  In  the  same  way,  if  some  additional 
feature,  such  as  the  ability  to  predict  nitrogen  oxides,  is  desired,  it  should  be  possible 
to  add  it  without  changing  the  core  of  the  remaining  mechanism. 

The  layout  of  the  approach  followed  in  this  work  is  provided  in  Fig.  15.  This  modu¬ 
lar  approach  relies  on  a  library  of  skeletal  mechanisms,  or  component  library,  to  build 
kinetic  schemes  for  mixtures  efficiently.  Given  the  type  of  fuel  and  the  application,  a  sur¬ 
rogate  composition  can  be  devised  through  the  optimization  process  presented  above.  On 
the  chemistry  side,  existing  detailed  mechanisms  for  individual  components  are  validated 
against  experimental  data,  then  reduced  for  various  conditions  to  the  skeletal  level  using 
two  reduction  techniques:  the  DRGEP  method,  and  a  chemical  lumping  technique.  The 
collection  of  skeletal  mechanisms  forms  the  so-called  component  library,  to  which  several 
specific  modules  can  be  added,  such  as  a  description  of  soot  formation  or  NOx  chemistry. 
The  surrogate  composition  and  its  future  use  dictate  the  choice  of  modules  that  have  to 
be  included  in  the  combined  skeletal  mechanism,  so  that  only  the  necessary  kinetics  are 
taken  into  account.  Then,  a  second  stage  of  reduction  is  applied  that  includes,  for  example, 
the  introduction  of  quasi-steady  state  assumptions.  The  obtained  reduced  model  can  be 
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Figure  15:  Layout  and  features  of  the  component  library  approach. 


validated  using  available  experimental  data.  The  results  of  this  last  step  will  provide  addi¬ 
tional  constraints  that  have  to  be  included  in  the  surrogate  definition  phase,  and  the  whole 
process  can  be  repeated  until  an  optimal  surrogate  is  obtained.  The  major  reduction  steps, 
that  is  DRGEP  and  lumping,  are  done  only  once,  which  minimizes  the  turn-over  time  to 
produce  the  surrogate  mechanism. 

2.3.3  Assumptions  and  Challenges 

The  component  library  approach,  and  the  efficiency  of  this  approach,  are  based  on  two 
major  assumptions.  The  first  one  is  that  two  different  large  hydrocarbon  molecules  will 
interact  during  combustion  only  at  the  level  of  small  radicals  and  decomposition  products 
that  are  present  already  in  both  detailed  mechanisms.  Cross-reactions  between  fuel-specific 
molecules  are  neglected.  The  main  argument  to  justify  this  assumption  involves  steric 
factors:  unless  the  pressure  is  very  high,  the  reactive  sites  of  large  molecules  or  radicals  are 
much  less  accessible  by  another  large  species  than  a  small  radical  such  as  OH.  The  validity 
of  this  assumption  was  shown  experimentally  by  Klotz  et  al  [34]  and  numerically  by  Zhao  et 
al  [35]  for  toluene/ alkane  mixtures.  In  case  the  cross- reactions  between  two  compounds  are 
shown  to  be  important,  incremental  sets  containing  the  relevant  reactions  could  be  added 


28 


in  the  library. 

The  second  assumption  is  that  not  all  reaction  pathways  are  important  for  all  config¬ 
urations.  This  assumption  allows  the  creation  of  modules  through  the  reduction  of  the 
detailed  mechanisms  on  complementary  sub-domains  of  the  parameter  space.  Then  smaller 
multi- component  mixtures  are  obtained  by  combining  only  the  modules  relevant  to  the  ap¬ 
plication.  The  best  example  of  this  segregated  approach  is  the  temperature  dependance  of 
the  hydrocarbon  oxidation  chemistry.  The  chemistry  can  be  separated  into  a  high  temper¬ 
ature  base  chemistry  and  a  low  temperature  module  that  can  be  superimposed  on  this  base 
chemistry  to  represent  the  characteristic  low  temperature  phenomena. 

The  chemistry  reduction  methods  and  how  to  determine  an  initial  surrogate  mixture 
composition  have  been  described  earlier.  How  the  various  modules  from  the  component 
library  and  incremental  sets  may  be  assembled  together  to  form  a  multi-component  mech¬ 
anism  remains  to  be  detailed.  This  step  is  simplified  considerably  by  the  fact  that  the 
detailed  mechanisms  have  been  reduced  already  to  a  skeletal  level;  however,  even  with 
this  simplification,  combining  kinetic  data  from  potentially  very  different  sources  is  a  non¬ 
trivial  task.  Species,  thermodynamic  and  transport  data,  and  elementary  reactions  need  to 
be  merged.  To  do  this  step  most  efficiently,  an  interactive  set-up  has  been  designed  that 
automatically  identifies  identical  species  and  reactions  from  the  modules  and  records  every 
multiple  choice  and  incompatibility  among  the  kinetic  data  sets.  For  example,  species  that 
have  the  same  formula  but  different  names  would  be  identified.  These  species  could  be 
identical  molecules,  different  isomers  of  the  same  species,  or  identical  reactions  with  sig¬ 
nificantly  different  rate  coefficients.  These  discrepancies  cannot  be  resolved  automatically; 
therefore,  they  are  subjected  to  the  user’s  expertise,  who  decides  which  option  is  the  best 
available. 

The  component  library  approach  described  above  has  been  used  in  several  applications. 
For  example,  the  modular  approach  is  illustrated  in  [36]  for  the  development  of  reduced 
models  for  a  gasoline  surrogate.  In  this  study,  specific  modules  for  low-temperature  ignition 
are  identified  and  added  to  a  reduced  model  describing  premixed  gasoline  fiames.  Blanquart 
et  al.  [37]  developed  a  kinetic  mechanism  describing  PAH  and  soot  formation,  to  which 
reduced  modules  for  heptane  and  iso-octane  oxidation  were  added  automatically.  A  reduced 
module  for  soot  formation  was  obtained  using  the  methods  described  above  and  used  to 
perform  direct  numerical  simulations  of  soot  formation  in  isotropic  turbulence  [38].  The 
next  section  illustrates  in  detail  the  development  of  a  reduced  model  for  a  jet  fuel  surrogate. 

2.3.4  Application  to  the  development  of  chemical  model  for  JP-8 

As  an  application  of  the  component  library  approach,  the  development  of  a  skeletal  model 
for  a  jet  fuel  surrogate  was  chosen.  Conventional  jet  fuels,  called  Jet-A  for  civil  applications 
and  JP-8  for  military,  have  a  relatively  high  average  carbon  number,  and  corresponding 
representative  hydrocarbons  are  difficult  to  handle  experimentally.  As  a  consequence,  little 
experimental  data  currently  exist  to  assist  in  the  development  of  detailed  chemical  mech¬ 
anisms  for  large  hydrocarbons.  As  new  data  become  available,  the  kinetic  models  will  be 
validated  and  refined  over  a  wider  range  of  conditions;  therefore,  the  present  work  does  not 
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attempt  to  produce  a  surrogate  model  that  will  match  perfectly  the  few  data  available  for 
jet  fuel.  Instead,  the  focus  is  placed  on  the  method  used  to  generate  the  reduced  model, 
especially  the  successive  validation  steps.  If  rigorous  validation  is  performed,  discrepancies 
in  the  results  can  be  attributed  mostly  to  inaccuracies  in  the  detailed  mechanisms  and  inad¬ 
equacy  of  the  initial  choice  of  components,  which  is  currently  very  limited.  The  subsequent 
development  of  better  surrogate  models  will  follow  from  improved  initial  experimental  and 
kinetic  data. 

Surrogate  Formulation  The  choice  of  the  individual  components  for  each  hydrocarbon 
class  was  motivated  mostly  by  the  availability  of  a  reasonably  validated  detailed  chemi¬ 
cal  mechanism.  They  include  7>heptane,  iso-octane,  and  n-dodecane  as  alkanes,  methyl- 
cyclohexane  as  naphthene,  and  benzene  and  toluene  as  aromatic  molecules.  The  base  mech¬ 
anism  is  taken  from  Blanquart  et  al.  37] .  This  mechanism  was  developed  specifically  to 
represent  soot  formation  and  was  validated  extensively  using  experimental  data  for  soot 
precursors,  PAHs,  and  soot  volume  fractions.  Mechanisms  for  7>heptane,  isooctane,  and 
methyl-cyclohexane  are  taken  from  Lawrence  Livermore  National  Laboratories  [39,  21,  40], 
and  n- dodecane  is  modeled  by  a  semi-detailed  mechanism  developed  by  Wang  et  al.  [41]. 

Three  surrogate  compositions  were  considered.  The  first  one  is  neat  7>dodecane.  The 
second  surrogate  was  developed  by  Violi  et  al.  [42,  43]  based  on  a  number  of  criteria,  includ¬ 
ing  sooting  tendency  and  distillation  characteristics.  The  third  one  is  obtained  through  the 
optimization  procedure  described  previously,  by  prescribing  the  hydrogen-to-carbon  ratio, 
the  average  carbon  number,  the  average  composition  in  terms  of  hydrocarbon  classes,  the 
cetane  number,  and  the  sooting  tendency  of  the  fuel.  Only  three  different  components  were 
included  in  the  optimization,  namely  n-dodecane,  methyl-cyclohexane,  and  toluene.  The 
resulting  composition  and  a  comparison  of  the  properties  of  the  three  surrogates  are  shown 
in  table  4. 


Average  Jet 
Fuel 

Neat 

Dodecane 

Surrog.  1 

Surrog.  2 

Composition 
[%  mol] 

Dodecane 

N/A 

100 

73.5 

45 

Iso-octane 

0 

5.5 

0 

MCH 

0 

10 

26.1 

Toluene 

0 

10 

28.9 

Benzene 

0 

1 

0 

H/C  ratio 

1.91 

2.17 

2.09 

1.09 

Formula 

Ci,H21 

Cl2H26 

C10  7H22.3 

C9.3H17.7 

Hydrocarbon 
composition 
[%  vol] 

Paraffins 

~60 

100 

88 

62 

Naphthenes 

~20 

0 

6.4 

20 

Aromatics 

~18 

0 

5.6 

18 

Cetane  number 

~42.7 

80 

73.4 

58 

Threshold  sooting  index 

~15 

5.2 

9.3 

16.3 

Table  4:  Compositions  for  jet  fuel  surrogate  used  in  comparison  with  experimental  measurements. 
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The  domain  of  applicability  of  the  reduced  surrogate  models  was  chosen  to  represent  the 
most  probable  operation  conditions  found  in  gas-turbine  engines,  namely  high  tempera  ture, 
atmospheric  to  high  pressures,  and  lean-to-rich  conditions.  Auto-ignition  and  premixed 
conditions  are  included  in  the  reduction  procedure,  as  newly  designed  engines  tend  to  run 
partially  premixed,  which  inevitably  leads  to  the  necessity  to  control  auto-ignition  at  oper¬ 
ating  conditions  [31].  The  detailed  mechanisms  for  the  individual  components  are  reduced 
independently  to  a  skeletal  level  using  DRGEP  and  chemical  lumping.  Then,  the  resulting 
modules  are  combined  together,  with  the  base  mechanism  of  Blanquart  et  al.  [37]  retained 
as  reference.  The  resulting  scheme  is  a  multi-component  mechanism  containing  only  181 
species,  able  to  accurately  describe  the  oxidation  of  n- heptane,  iso-octane,  ?>dodecane, 
methyl-cyclohexane,  benzene,  toluene,  and  the  soot  precursor  acetylene.  Validation  is  per¬ 
formed  at  each  stage,  and  includes: 

•  Comparison  of  the  detailed  mechanism  with  experimental  data 

•  Comparison  of  simulations  obtained  using  the  single  component  skeletal  mechanism 
and  the  detailed  mechanism  from  which  it  is  extracted 

•  Comparison  of  simulations  obtained  for  each  individual  component  using  the  combined 
skeletal  mechanism  and  the  corresponding  single  component  detailed  mechanism 

•  Assessment  of  the  performance  of  the  combined  surrogate  mechanism  by  comparing 
simulation  results  with  experimental  data  for  jet  fuels. 

Comparisons  with  experimental  data  are  described  in  detail  in  [20].  The  pure  com¬ 
ponents  mentioned  above  were  extensively  studied  in  a  wide  range  of  configurations,  such 
as  auto-ignition,  plug-flow  reactors,  premixed  and  counter-flow  flames  for  various  tempera¬ 
tures,  pressure  and  equivalence  ratios.  The  multi-component  model  showed  good  agreement 
with  both  the  original  detailed  schemes  and  the  various  experimental  sets. 

Simulation  results  for  ignition  delay  times  of  jet  fuel  are  shown  in  Fig.  16.  Figure  16(a) 
compares  two  different  pressures,  20  and  50  atm,  Fig.  16(b)  compares  two  levels  of  dilution 
for  stoichiometric  mixtures  of  jet  fuel  and  air,  and  Fig.  16(c)  compares  ignition  delay  times 
for  two  different  equivalence  ratios,  0.5  and  1.  The  conclusions  are  similar  in  each  of 
the  cases.  First,  all  surrogates  over-predict  auto-ignition  by  a  factor  two  to  three.  New 
experimental  data  from  Vasu  et  al.  [44]  show  that  7>dodecane/air  mixtures  ignition  delay 
times  are  close  to  those  of  jet  fuel.  As  7>dodecane  is  the  major  component  in  all  surrogates, 
this  observation  tends  to  indicate  that  the  7>dodecane  mechanism  does  not  predict  auto- 
ignition  well.  Also,  rc-dodecane  is  the  component  that  ignites  the  fastest.  As  surrogate  #2 
contains  less  of  the  large  alkane  than  the  others,  this  surrogate  shows  the  biggest  discrepancy 
compared  to  experiments  at  high  temperature. 

A  second  observation  concerns  the  onset  of  the  NTC  region.  Experimental  data  show 
no  NTC  behavior  for  temperatures  larger  than  850  K.  However,  both  pure  7>dodecane 
and  Surrogate  #1  have  a  strong  NTC  behavior  for  temperatures  around  1000  K.  This 
NTC  is  much  weaker  for  Surrogate  #2.  Disregarding  the  shift  due  to  the  slow  ignition  of 
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(c)  Effect  of  equivalence  ratio  (20.7%  O2,  20 
atm). 


Figure  16:  Auto-ignition  of  JetA/JP-8.  Comparison  between  experimental  data  (squares,  [44]),  and  simula¬ 
tion  results  obtained  with  pure  dodecane  (dash-dotted  lines),  surrogate  #1  (dashed  lines)  and  surrogate  #2 
(solid  lines). 


r^dodecane  in  the  simulation,  Surrogate  #2  better  represents  the  jet  fuel  experimental  data. 

Experimental  and  simulated  species  concentration  profiles  obtained  in  a  premixed  at¬ 
mospheric  kerosene  Hame,  with  an  initial  temperature  of  473  K  [45],  are  shown  in  Fig.  17. 
Prediction  of  oxygen  consumption  and  major  product  formation,  such  as  CO  and  CO2,  is  in 
acceptable  agreement  with  the  measurements.  The  precise  compositions  of  the  surrogates 
have  very  little  effect  on  the  evolution  of  the  major  species;  however,  the  concentration  of 
specific  intermediate  species  clearly  are  impacted  by  the  initial  composition.  For  exam¬ 
ple,  neat  7vdodecane  tends  to  produce  very  little  benzene,  whereas  the  presence  of  toluene 
and  benzene  from  the  inlet  increases  drastically  the  maximum  concentration  of  benzene  in 
the  flame.  Other  common  intermediates  that  are  natural  decomposition  products  of  any 
sufficiently  large  hydrocarbon  fuel  are  produced  in  similar  quantity  by  the  three  surrogates. 
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Figure  17:  Premixed  burner-stabilized  kerosene  flame.  Comparison  between  experimental  data 
(squares,  [45]),  and  simulation  results  obtained  with  puredodecane  (dash-dotted  lines),  surrogate  #1  (dashed 
lines)  and  surrogate  #2  (solid  lines). 

The  present  experimental  validation  set  is  not  sufficient  to  conclude  with  certitude  on 
the  performances  of  the  tested  surrogates  except,  perhaps,  for  the  neat  n-dodecane  case 
that  is  clearly  not  representative  of  either  jet  fuel  auto-ignition  or  flame  oxidation.  Indeed, 
too  much  uncertainty  still  exists  in  the  detailed  kinetic  models  considered  in  this  study; 
however,  should  any  of  these  mechanisms  be  improved  in  the  near  future,  only  a  few  steps 
need  to  be  repeated,  namely  the  reduction  of  this  particular  mechanism,  combination  with 
the  other  skeletal  schemes,  and  validation  steps  for  all  components.  Provided  that  the 
reference  mechanism  is  left  unchanged,  this  latter  stage  should  not  introduce  any  problems 
in  the  untouched  components,  and  the  whole  process  is  expected  to  be  quite  fast. 
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3  Conclusions 


In  this  project,  we  have  focussed  on  two  topics:  advanced  combustion  models  for  large* 
eddy  simulations  (LES)  and  reduced  chemical  mechanisms  for  JP-8  surrogate  fuels.  We  have 
made  substantial  progress  on  both  of  these  topics.  The  aim  of  the  combustion  LES  modeling 
part  was  to  advance  the  models  for  non-premixed  and  premixed  combustion  towards  a 
generalized  combustion  model  that  covers  all  combustion  regimes.  Towards  this  end,  we 
have  made  three  major  model  advancements.  First,  a  dynamic  model  for  the  turbulent 
burning  velocity  was  developed.  This  model  eliminates  the  use  of  model  constants  in  the 
G-equation  formalism,  which  is  used  to  compute  the  flame  position.  The  propagation  term 
appearing  in  this  equation  now  is  computed  dynamically.  The  model  formulation  comes 
from  applying  the  constraint  that  the  model  must  produce  consistently  the  same  flame 
behavior  at  all  LES  filter  levels.  Additionally,  a  new  filtering  approach  has  been  developed 
that  now  allows  for  the  use  of  the  volumetric  filters.  This  approach  enables  explicit  filtering 
of  the  flame  on  the  LES  test  filter  length  scale.  Direct  numerical  simulations  (DNS)  of 
flame  propagation  in  a  turbulent  flow  were  performed  to  assess  the  model.  The  model  was 
found  to  predict  the  DNS  results  very  accurately.  Advantages  of  the  new  model  were  found 
specifically  in  the  early  development  of  the  turbulent  flame  before  a  statistical  steady  state 
is  reached. 

Second,  a  flame  structure  model  was  developed.  This  additional  model  is  necessary 
because  the  solution  of  the  G-equation  provides  only  the  mean  position  of  a  given  tempera¬ 
ture  iso-surface,  typically  chosen  to  be  that  of  the  inner  layer  temperature.  For  a  sub-filter 
Damkohler  number  larger  than  unity,  the  entire  flame  appears  on  the  sub-filter  scale  and 
requires  no  further  treatment;  however,  for  small  Damkohler  numbers,  the  flame  is  resolved 
only  partially  and  cannot  be  described  by  the  G-equation  alone.  In  the  flame  structure 
model,  an  equation  for  the  reaction  progress  variable  is  solved  in  addition  to  the  G-equation. 
This  equation  shows  a  sharp  jump  across  the  filtered  flame  and  therefore  cannot  be  used  to 
track  the  flame  accurately.  In  the  new  formulation,  both  the  transport  term  and  the  chemi¬ 
cal  source  term  are  modeled  such  that  the  filtered  position  of  the  inner  layer  predicted  from 
the  progress  variable  equation  is  consistent  with  that  of  the  G-equation.  This  treatment 
ensures  the  correct  propagation  speed,  and  the  filtered  temperature  and  density  then  can  be 
computed  from  the  progress  variable  field.  Our  models  for  non-premixed  combustion  also 
use  a  reaction  progress  variable  equation.  The  progress  variable  therefore  provides  a  way 
to  merge  the  two  models  for  premixed  and  non-premixed  combustion.  The  source  terms, 
however,  must  be  modeled  as  a  function  of  the  combustion  regime. 

Third,  a  formalism  to  detect  the  combustion  regime  locally  was  developed.  The  so- 
called  flame  index  was  proposed  previously,  but  we  have  shown  that  this  quantity  is  too 
simplistic  and  does  not  work  in  general.  Our  method  is  based  on  an  asymptotic  analysis 
of  the  limiting  cases.  The  existence  of  the  flame  structure  model  and  a  combustion  regime 
identifier  provide  a  basis  for  developing  a  consistent  generalized  model  that  works  across 
different  combustion  regimes. 

In  the  chemistry  part  of  the  project,  we  have  focussed  on  two  issues,  the  component  li- 
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brary  approach  and  the  automatic  reduction  procedure.  In  the  component  library  approach, 
reduced  chemical  mechanisms  are  stored  as  modules  in  a  component  library.  Several  mod¬ 
ules  could  be  included  in  the  component  library  for  each  fuel  component,  such  as  a  basis 
module  for  high-temperature  chemistry  and  additional  incremental  sets  for  low-temperature 
chemistry  and  soot  precursors.  In  addition  to  populating  the  component  library  with  sev¬ 
eral  components,  we  have  made  the  reduction  process  more  automatic  and  have  developed 
tools  for  the  automatic  combination  of  mechanisms  comprised  of  different  components.  The 
automatic  reduction  procedure  itself  was  developed  by  refining  the  DRGEP  method.  This 
method  now  is  used  to  eliminate  species  and  reactions.  Additionally,  we  have  introduced 
an  automatic  procedure  to  select  steady  state  species,  and  we  have  developed  a  new  ac¬ 
curate  method  for  species  lumping.  The  resulting  automatic  reduction  procedure  is  based 
on  the  combination  of  all  these  approaches.  The  method  was  tested  on  the  largest  avail¬ 
able  mechanisms  for  various  hydrocarbon  molecules,  such  as  n-heptane,  iso-octane,  and 
methylcyclohexane.  The  mechanisms  were  reduced  automatically,  and  the  component  li¬ 
brary  approach  was  validated  successfully.  The  resulting  mechanism  is  quite  small  compared 
to  the  original,  but  similarly  small  mechanisms  have  been  achieved  using  other  reduction 
methods.  What  is  unique  about  our  approach  is  that,  in  contrast  to  other  methods,  the 
entire  reduction  procedure  is  automatic  and  based  on  a  single  evaluation  of  the  production 
rates  with  the  detailed  mechanism.  The  reduction,  therefore,  is  cheap  and  occurs  all  at  once. 
This  feature  is  important,  since  it  could  form  the  basis  of  an  automatic  adaptive  chemistry 
approach,  where  the  reduction  is  performed  automatically  during  a  reactive  computational 
fluid  dynamics  simulation. 
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